rm(list=ls())
library(data.table)
library(ggplot2)
library(polywog)
library(MASS)
library(doParallel)
library(directlabels)
load("ensemble.Rdata")
ensemble$sddistance <- ensemble$sd.distance/.12
ensemble <- ensemble[, sd.distance:=NULL]
ensemble <- ensemble[, diff:=parties-mean.num.gov.parties]
ensemble <- ensemble[,alpha5:=ifelse(alpha>.5, 1, 0)]

ensemble$sd.eccentricity <- ensemble$se.eccentricity*sqrt(500)
ensemble$sd.representation <- ensemble$se.representation*sqrt(500)
ensemble$sd.gov.representation2 <- ensemble$se.gov.representation2*sqrt(500)


#########################
# FUNCTION FOR PLOTTING #
#########################
plot.vars <- function (data, varying, by="",
                         parties.val=5,
                         discount.val=.5,
                         alpha.val=.5,
                         sd.distance.val=1
) {
 if (varying=="parties"|by=="parties") {
     parties.val <- unique(data[,parties])
 }
 if (varying=="discount"|by=="discount") {
     discount.val <- unique(data[,discount])
  }
 if (varying=="sd.distance"|by=="sd.distance") {
     sd.distance.val <- unique(data[,sd.distance])
 }
 if (varying=="rule"|by=="rule") {
     rule.val <- unique(data[,rule])
 }
 if (varying=="alpha"|by=="alpha") {
     alpha.val <- unique(data[,alpha])
 }
 output <- data[round==76&
                   discount %in% discount.val
                 & parties %in% parties.val
                 & alpha %in% alpha.val
                 & sd.distance %in% sd.distance.val]
 output <- as.data.frame(output)
 plot.data <-data.frame(x=output[,varying],
                          y=output$mean,
                          lower=output$mean-(1.96*output$se),
                          upper=output$mean+(1.96*output$se),
                          variable=output$variable,
                          rule=output$rule)
 if (by!="") {
   plot.data$z <- output[,by]
 } 
     p <- ggplot(plot.data, aes(x=x, y=y, group=rule))
     #p <- p + geom_segment(aes(colour=variable, y=upper, yend=lower, xend=x))
     #p <- p + geom_point(aes(colour=variable))
     p <- p + geom_ribbon(fill="grey73", aes(ymin=lower, ymax=upper))
     p <- p + geom_line(aes(colour=as.factor(rule)))
     p <- p + theme_bw()
     p <- p + theme(panel.grid.major.x = element_blank())
     p <- p + theme(panel.grid.major.y = element_blank())
     p <- p + theme(panel.grid.minor.x = element_blank())
     p <- p + theme(panel.grid.minor.y = element_blank())
     p <- p + xlab(varying)
     p <- p + ylab("")
     p <- p + scale_colour_manual("Rules",
                                  breaks=c(2, 3, 4, 7, 8),
                                  values=c("red", "blue", "yellow", "green", "pink"),
                                  labels=c("Aggregator", "Sat Governator", "Sticker",
                                        "Hunter", "Hunt Governator"))
     ## p <- p + scale_colour_manual("Rules", labels=c("Aggregator", "Sat Governator", "Sticker",
     ##                                    "Hunter", "Hunt Governator"))
     if (by!="") {
         p <- p + facet_grid(variable~z, scales="free")
      } else {
         p <- p + facet_grid(variable~., scales="free")
      }
return(p)
} #end plot.vars

#This function data-mines using polywog but eliminates terms from the equation
polywog.select <- function(formula, data, degrees.cv, boot=5, select=.05, lambda=NULL) {
    cv <- cv.polywog(formula=as.formula(formula),
                     data=as.data.frame(data),
                     degrees.cv=as.numeric(degrees.cv),
                     .parallel=TRUE)
    fit <- cv$polywog.fit
    fit <- bootPolywog(fit, nboot = boot, .parallel=TRUE)  
    sig <- which((abs(summary(fit)$coefficients[,1]/summary(fit)$coefficients[,2])>1.96)==TRUE)
    rhs <- names(sig)
    rhs <- gsub("[.]", ")*I(", rhs) 
    if (rhs[1]=="(Intercept)") {
        rhs <- rhs[-1]
    } else {
        rhs[1] <- "-1"
    }
    if (length(rhs)>1) {
        rhs <- ifelse(rhs!="-1", paste("I(", rhs, ")", sep=""), "-1")
        rhs <- paste(rhs, collapse="+")
    } else if (length(rhs)==1) {
        rhs <- paste("I(", rhs, ")", sep="")
    } else if (length(rhs)==0){
      rhs <- "1"
    }  
    new.formula <- as.formula(paste(as.character(formula[2]), "~", rhs))
    return(lm(new.formula, data))
}

#This function makes predictions based on polywog.select results.
predict.polywog <- function(pmodel, parties=3, discount=0, sddistance=0, alpha=0, rule="") {
 input.data <- expand.grid(0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                           parties, discount, sddistance, alpha)
 names(input.data) <- c("sd.eccentricity", "sd.representation", "sd.gov.representation2",
                        "mean.eccentricity", "mean.representation", "mean.gov.representation",
                        "mean.gov.representation2",
                        "mean.num.gov.parties", "mean.gov.seats", "mean.gov.cohesion",
                        "mean.oversized", "mean.majority", "mean.ENgov", 
                        "parties", "discount", "sddistance", "alpha")
 input.data <- as.data.frame(input.data)
 output.data <- as.data.frame(input.data)
 input.data$alpha1 <- ifelse(input.data$alpha==1, 1, 0)
 input.data$alpha0 <- ifelse(input.data$alpha==0, 1, 0)
 input.data$alpha5 <- ifelse(input.data$alpha==.5, 1, 0)
 input.data$discount1 <- ifelse(input.data$discount==1, 1, 0)
 input.data$discount0 <- ifelse(input.data$discount==0, 1, 0)
 input.data$sddistance0 <- ifelse(input.data$sddistance==0, 1, 0)                         
 input <- model.matrix(pmodel, data=input.data)
 betas <- mvrnorm(1000, coef(pmodel), vcov(pmodel))
 values <- input %*% t(betas)
 means <- apply(values, 1, mean)
 lower <- apply(values, 1, quantile, p=.025)
 upper <- apply(values, 1, quantile, p=.975)
 return(as.data.table(cbind(output.data, lower, means, upper,
                            type=as.character(pmodel$terms[[2]]),
                            rule=rule)))
}



#This function plots based on predict polywog.
plot.polywog <- function(x, varying, by="", plot=TRUE, color=FALSE, interest="eccentricity",
                         original=NULL) {
 interest.mean <- paste0("mean.", interest)
 interest.se <- paste0("se.", interest)
 x <- as.data.table(x)
 if (plot==TRUE) {
   x <- as.data.frame(x)
     if (by!="") {
       plot.data <- data.frame(x=x[,varying],
                               y=x$means,
                               lower=x$lower,
                               upper=x$upper,
                               z=x[,by],
                               original=0)
#      plot.data <- as.data.frame(cbind(x,y, y.se, z))
      } else {    
       plot.data <- data.frame(x=x[,varying],
                               y=x$means,
                               lower=x$lower,
                               upper=x$upper,
                               original=0)
   }
   if (is.null(original)==FALSE) {
       #original <- as.data.frame(original)
      if (varying=="parties") {
          parties.val <- unique(c(original[,parties], x[,"parties"]))
      } else {
          parties.val <- unique(x[,"parties"])
      }
      if (varying=="discount") {
          discount.val <- unique(c(original[,discount], x[,"discount"]))
      } else {
          discount.val <- unique(x[,"discount"])
      }
      ## if (varying==memory) {
      ##     memory.val <- unique(c(original[,memory], x[,memory]))
      ## } else {
      ##     memory.val <- unique(x[,memory])
      ## }
      ## if (varying==evol) {
      ##     evol.val <- unique(c(original[,evol], x[,evol]))
      ## } else {
      ##     evol.val <- unique(x[,evol])
      ## }
      if (varying=="sddistance") {
          sddistance.val <- unique(c(original[,sddistance], x[,"sddistance"]))
      } else {
          sddistance.val <- unique(c(x[,"sddistance"], 0.4175))
      }
      if (varying=="alpha") {
          alpha.val <- unique(c(original[,alpha], x[,"alpha"]))
      } else {
          alpha.val <- unique(x[,"alpha"])
      }
 output <- original[discount %in% discount.val& parties %in%parties.val  & alpha %in% alpha.val & sddistance %in% sddistance.val]
       output <- as.data.frame(output)
       if (by!="") {
        plot.data <- rbind(plot.data, 
                                cbind(x=output[,varying],
                                y=output[,interest.mean],
                               lower=output[,interest.mean]-(1.96*output[,interest.se]),
                               upper=output[,interest.mean]+(1.96*output[,interest.se]),
                                z=output[,by],
                                original=1)
                           )
#       plot.data <- as.data.frame(cbind(x,y, y.se, z))
       } else {    
        plot.data <- rbind(plot.data, 
                                cbind(x=output[,varying],
                                y=output[,interest.mean],
                               lower=output[,interest.mean]-(1.96*output[,interest.se]),
                               upper=output[,interest.mean]+(1.96*output[,interest.se]),
                                original=1)
                           )
      }
     }
     p <- ggplot(plot.data)
     p <- p + theme_bw()
     p <- p + theme(panel.grid.major.x = element_blank())
     p <- p + theme(panel.grid.major.y = element_blank())
     p <- p + theme(panel.grid.minor.x = element_blank())
     p <- p + theme(panel.grid.minor.y = element_blank())
     p <- p + xlab(varying)
     p <- p + ylab(interest.mean)
     if (by!="") {
      if (color==FALSE) {
       p <- p+ facet_grid(~z)
       p <- p + geom_ribbon(aes(x=x,
                                ymin=lower,
                                ymax=upper),
                            fill="grey74",
                            data=plot.data[plot.data$original==0,])
       p <- p + geom_line(aes(x=x, y=y),
                            data=plot.data[plot.data$original==0,])
      } else {
       p <- p + geom_ribbon(aes(x=x,
                                ymin=lower,
                                ymax=upper),
                            fill="grey74",
                            data=plot.data[plot.data$original==0,])
       p <- p + geom_line(aes(x=x, y=y, color=as.factor(z)),
                            data=plot.data[plot.data$original==0,])
      }
     } else {
         p <- p + geom_ribbon(aes(x=x,
                                 ymin=lower,
                                 ymax=upper),
                             fill="grey74",
                            data=plot.data[plot.data$original==0,])
         p <- p + geom_line(aes(x=x, y=y),
                            data=plot.data[plot.data$original==0,])
     }
     if (is.null(original)==FALSE) {
        p <- p + geom_segment(aes(x = x,
                                  y = lower,
                                  xend = x,
                                  yend = upper),
                              size=.1, data=plot.data[plot.data$original==1,])
       p <- p + geom_point(aes(x = x, y = y),
                           data=plot.data[plot.data$original==1,])
       p <- p+ facet_grid(~z)
     }
     return(p)
 } else {
     return(output)
 }
}


####################
# PLOT REAL VALUES #
####################

## parties:       3    4    5    6  7  8
## aspiration: 0.25 0.50 0.75 0.90
## memory:      2  3  4  5  6  7  8  9 10 
## sd.distance: 0.00 0.25 0.50 1.00 1.50 2.00
## discount:    0.00 0.10 0.25 0.50 0.75 0.90 1.00
## alpha:       0.00 0.10 0.25 0.50 0.75 0.90 1.00 (1=policy-seeking)

## plot.vars(ensemble, "parties", 
##           parties=5, discount=.5, alpha=.5, sd.distance=1)
## plot.vars(ensemble, "sd.distance", by="alpha", parties=5, discount=.5)
## plot.vars(ensemble, "discount","alpha", parties=5, discount=.5)
## plot.vars(ensemble, "alpha", parties=5, discount=.5, sd.distance=1)
## plot.vars(ensemble, "parties","alpha",
##           parties=5, discount=.5, alpha=.5, sd.distance=1)
## plot.vars(ensemble, "parties","sd.distance",
##           parties=5, discount=.5, alpha=.5, sd.distance=1)
## plot.vars(ensemble, "alpha", "sd.distance")



###########################################
# POLYWOG RESULTS, PREDICT, AND PLOT THEM #
###########################################


## cl <- makeCluster(4)
## registerDoParallel(cl)

## ##############
## # Aggregator #
## ##############
## print("Aggregator")
## set.seed(83758) #from random.org
## data=ensemble[ensemble$rule==2&ensemble$round==76,]
## formula=as.formula(sd.eccentricity~
##                       alpha
##                       +parties
##                       +discount
##                       +sddistance)
## A.ecc <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=40)
## save(A.ecc, file="sdA.ecc.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(sd.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## A.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(A.rep, file="sdA.rep.Rdata")

## set.seed(83758) #from random.org
## formula=as.formula(sd.gov.representation2~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## A.gov.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(A.gov.rep, file="sdA.gov.rep2.Rdata")

## ##################
## # Sat Governator #
## ##################
## print("SatGov")
## set.seed(83758) #from random.org
## data=ensemble[ensemble$rule==3&ensemble$round==76,]
## formula=as.formula(sd.eccentricity2~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.ecc <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.ecc, file="sdSG.ecc.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(sd.representation2~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.rep, file="sdSG.rep.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(sd.gov.representation2~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## SG.gov.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(SG.gov.rep, file="sdSG.gov.rep2.Rdata")

## ###########
## # Sticker #
## ###########
## print("Sticker")
## set.seed(83758) #from random.org
## data=ensemble[ensemble$rule==4&ensemble$round==76,]
## formula=as.formula(sd.eccentricity~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.ecc <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.ecc, file="sdS.ecc.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(sd.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.rep, file="sdS.rep.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(sd.gov.representation2~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## S.gov.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(S.gov.rep, file="sdS.gov.rep2.Rdata")


## ## HUNTER #
## print("Hunter")
## set.seed(83758) #from random.org
## data=ensemble[ensemble$rule==7]
## formula=as.formula(sd.eccentricity~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.ecc <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.ecc, file="sdH.ecc.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(sd.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.rep, file="sdH.rep.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(sd.gov.representation2~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## H.gov.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(H.gov.rep, file="sdH.gov.rep.Rdata")


## ###################
## # Hunt Governator #
## ###################
## print("HunGov")
## set.seed(83758) #from random.org
## data=ensemble[ensemble$rule==8&ensemble$round==76,]
## formula=as.formula(sd.eccentricity~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## HG.ecc <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(HG.ecc, file="sdHG.ecc.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(sd.representation~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## HG.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(HG.rep, file="sdHG.rep.Rdata")


## set.seed(83758) #from random.org
## formula=as.formula(sd.gov.representation2~
##                       alpha+alpha5
##                       +parties
##                       +discount
##                       +sddistance)
## HG.gov.rep <- polywog.select(formula=formula, data=data, degrees.cv=1:5, boot=400)
## save(HG.gov.rep, file="sdHG.gov.rep.Rdata")
## stopCluster(cl)



load("sdA.ecc.Rdata")
load("sdA.rep.Rdata")
load("sdA.gov.rep2.Rdata")
load("sdSG.ecc.Rdata")
load("sdSG.rep.Rdata")
load("sdSG.gov.rep2.Rdata")
load("sdS.ecc.Rdata")
load("sdS.rep.Rdata")
load("sdS.gov.rep2.Rdata")
load("sdH.ecc.Rdata")
load("sdH.rep.Rdata")
load("sdH.gov.rep.Rdata")
load("sdHG.ecc.Rdata")
load("sdHG.rep.Rdata")
load("sdHG.gov.rep2.Rdata")

v.parties <- seq(3,8,1)
v.discount <- .5
v.sddistance <- 1
v.alpha <- .5
Ae <- predict.polywog(A.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Se <- predict.polywog(S.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGe <- predict.polywog(HG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGe <- predict.polywog(SG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
He <- predict.polywog(H.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ar <- predict.polywog(A.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sr <- predict.polywog(S.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGr <- predict.polywog(HG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGr <- predict.polywog(SG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hr <- predict.polywog(H.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ag <- predict.polywog(A.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sg <- predict.polywog(S.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
SGg <- predict.polywog(SG.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hg <- predict.polywog(H.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
HGg <- predict.polywog(HG.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
t.data <- data.frame(x=c(7.5, 7.5, 7.6, 6.9, 4),
                     y=c(4.35, 3.85, 2.10, 1.17, .95),
                     label=c("S", "A", "H",
                         "SG", "Governator"),
                     type="eccentricity.mean")
data <- rbindlist(list(Ae, Se, He, SGe, HGe, Ar, Sr, HGr, SGr, Hr, Ag, Sg, Hg, SGg, HGg))
data$type <- as.character(data$type)
data$type <- ifelse(data$type=="sd.eccentricity", "Eccentricity", data$type)
data$type <- ifelse(data$type=="sd.representation", "System Misery", data$type)
data$type <- ifelse(data$type=="sd.gov.representation", "Government Misery", data$type)
data$type <- ifelse(data$type=="sd.gov.representation2", "Government Misery", data$type)
data$type <- factor(data$type, levels=c("Eccentricity", "System Misery", "Government Misery"))
data <- data[,label.x:=max(parties)+.2, list(rule, type)]
data <- data[,label.y:=means[parties+.2==label.x], list(rule, type)]
data$label.y <- ifelse(data$type=="Eccentricity"&data$rule=="G", data$label.y-.001,data$label.y)
data$label.y <- ifelse(data$type=="Eccentricity"&data$rule=="S", data$label.y+.001,data$label.y)
data$label.y <- ifelse(data$type=="System Misery"&data$rule=="S", data$label.y-.0002,data$label.y)
data$label.y <- ifelse(data$type=="System Misery"&data$rule=="H", data$label.y-.0002,data$label.y)
p <- ggplot(data=data, aes(x=parties, y=means, group=rule))
p <- p + geom_ribbon(aes(ymin=lower, ymax=upper, x=parties), fill="gray80")
p <- p + facet_grid(type~., scales="free")
p <- p + geom_line(aes(colour=as.factor(rule)))
p <- p + theme_bw()
p <- p + theme(panel.grid.major.x = element_blank())
p <- p + theme(panel.grid.major.y = element_blank())
p <- p + theme(panel.grid.minor.x = element_blank())
p <- p + theme(panel.grid.minor.y = element_blank())
p <- p + ylab("Expected Mean of Variability")
p <- p + xlab("Number of Parties")
p <- p + geom_text(aes(x=label.x, y=label.y, label=rule))
## p <- p + scale_colour_manual(name="Rule", values=c("red", "gold2", "green", "blue", "magenta"))
p <- p + coord_fixed()
## p <- p + geom_dl(aes(label=as.factor(rule)), method = list(dl.trans(x = x + 0.2), "last.points", cex = 0.8))
p <- p + theme(legend.position="none")
p <- p + theme(text = element_text(size=15))
p
ggsave("sd_parties.png", dpi=1000)


v.parties <- 5
v.discount <- .5
v.sddistance <- 1
v.alpha <- seq(0,1,.1)
Ae <- predict.polywog(A.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Se <- predict.polywog(S.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGe <- predict.polywog(HG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGe <- predict.polywog(SG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
He <- predict.polywog(H.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ar <- predict.polywog(A.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sr <- predict.polywog(S.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGr <- predict.polywog(HG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGr <- predict.polywog(SG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hr <- predict.polywog(H.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ag <- predict.polywog(A.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sg <- predict.polywog(S.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
SGg <- predict.polywog(SG.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hg <- predict.polywog(H.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
HGg <- predict.polywog(HG.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
t.data <- data.frame(x=c(7.5, 7.5, 7.6, 6.9, 4),
                     y=c(4.35, 3.85, 2.10, 1.17, .95),
                     label=c("S", "A", "H",
                         "SG", "Governator"),
                     type="eccentricity.mean")
data <- rbindlist(list(Ae, Se, He, SGe, HGe, Ar, Sr, HGr, SGr, Hr, Ag, Sg, Hg, SGg, HGg))
data$type <- as.character(data$type)
data$type <- ifelse(data$type=="sd.eccentricity", "Eccentricity", data$type)
data$type <- ifelse(data$type=="sd.representation", "System Misery", data$type)
data$type <- ifelse(data$type=="sd.gov.representation", "Government Misery", data$type)
data$type <- ifelse(data$type=="sd.gov.representation2", "Government Misery", data$type)
data$type <- factor(data$type, levels=c("Eccentricity", "System Misery", "Government Misery"))
data <- data[,label.x:=max(alpha)+.02, list(rule, type)]
data <- data[,label.y:=means[alpha+.02==label.x], list(rule, type)]
data$label.y <- ifelse(data$type=="Eccentricity"&data$rule=="S", data$label.y+.001,data$label.y)
data$label.y <- ifelse(data$type=="System Misery"&data$rule=="G", data$label.y+.0002,data$label.y)
data$label.y <- ifelse(data$type=="Government Misery"&data$rule=="G",
                       data$label.y+.004,data$label.y)
p <- ggplot(data=data, aes(x=alpha, y=means, group=rule))
p <- p + geom_ribbon(aes(ymin=lower, ymax=upper, x=alpha), fill="gray80")
p <- p + facet_grid(type~., scales="free")
p <- p + geom_line(aes(colour=as.factor(rule)))
p <- p + theme_bw()
p <- p + theme(panel.grid.major.x = element_blank())
p <- p + theme(panel.grid.major.y = element_blank())
p <- p + theme(panel.grid.minor.x = element_blank())
p <- p + theme(panel.grid.minor.y = element_blank())
p <- p + ylab("Expected Mean of Variability")
p <- p + xlab("Policy-Seekingness")
p <- p + geom_text(aes(x=label.x, y=label.y, label=rule))
## p <- p + scale_colour_manual(name="Rule", values=c("red", "gold2", "green", "blue", "magenta"))
p <- p + coord_fixed()
## p <- p + geom_dl(aes(label=as.factor(rule)), method = list(dl.trans(x = x + 0.2), "last.points", cex = 0.8))
p <- p + theme(legend.position="none")
p <- p + theme(text = element_text(size=15))
p
ggsave("sd_alpha.png", dpi=1000)


v.parties <- 5
v.discount <- .5
v.sddistance <- seq(0,2, .1)
v.alpha <- .5
Ae <- predict.polywog(A.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Se <- predict.polywog(S.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGe <- predict.polywog(HG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGe <- predict.polywog(SG.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
He <- predict.polywog(H.ecc, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ar <- predict.polywog(A.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sr <- predict.polywog(S.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
HGr <- predict.polywog(HG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
SGr <- predict.polywog(SG.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hr <- predict.polywog(H.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
Ag <- predict.polywog(A.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="A")
Sg <- predict.polywog(S.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="S")
SGg <- predict.polywog(SG.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="SG")
Hg <- predict.polywog(H.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="H")
HGg <- predict.polywog(HG.gov.rep, parties=v.parties, discount=v.discount, sddistance=v.sddistance,
                     alpha=v.alpha, rule="G")
t.data <- data.frame(x=c(7.5, 7.5, 7.6, 6.9, 4),
                     y=c(4.35, 3.85, 2.10, 1.17, .95),
                     label=c("S", "A", "H",
                         "SG", "Governator"),
                     type="eccentricity.mean")
data <- rbindlist(list(Ae, Se, He, SGe, HGe, Ar, Sr, HGr, SGr, Hr, Ag, Sg, Hg, SGg, HGg))
data$type <- as.character(data$type)
data$type <- ifelse(data$type=="sd.eccentricity", "Eccentricity", data$type)
data$type <- ifelse(data$type=="sd.representation", "System Misery", data$type)
data$type <- ifelse(data$type=="sd.gov.representation", "Government Misery", data$type)
data$type <- ifelse(data$type=="sd.gov.representation2", "Government Misery", data$type)
data$type <- factor(data$type, levels=c("Eccentricity", "System Misery", "Government Misery"))
data <- data[,label.x:=max(sddistance)+.03, list(rule, type)]
data <- data[,label.y:=means[sddistance+.03==label.x], list(rule, type)]
## data$label.y <- ifelse(data$type=="Eccentricity"&data$rule=="H", data$label.y+.002,data$label.y)
## data$label.y <- ifelse(data$type=="System Misery"&data$rule=="S", data$label.y+.0003,data$label.y)
## data$label.y <- ifelse(data$type=="Government Misery"&data$rule=="G",
##                        data$label.y+.005,data$label.y)
## data$label.y <- ifelse(data$type=="Government Misery"&data$rule=="SG",
##                        data$label.y-.005,data$label.y)
## data$label.y <- ifelse(data$type=="Government Misery"&data$rule=="A",
##                        data$label.y-.002,data$label.y)
p <- ggplot(data=data, aes(x=sddistance, y=means, group=rule))
p <- p + geom_ribbon(aes(ymin=lower, ymax=upper, x=sddistance), fill="gray80")
p <- p + facet_grid(type~., scales="free")
p <- p + geom_line(aes(colour=as.factor(rule)))
p <- p + theme_bw()
p <- p + theme(panel.grid.major.x = element_blank())
p <- p + theme(panel.grid.major.y = element_blank())
p <- p + theme(panel.grid.minor.x = element_blank())
p <- p + theme(panel.grid.minor.y = element_blank())
p <- p + ylab("Expected Mean of Variability")
p <- p + xlab("Ideal Point Dispersion")
p <- p + geom_text(aes(x=label.x, y=label.y, label=rule))
## p <- p + scale_colour_manual(name="Rule", values=c("red", "gold2", "green", "blue", "magenta"))
p <- p + coord_fixed()
## p <- p + geom_dl(aes(label=as.factor(rule)), method = list(dl.trans(x = x + 0.2), "last.points", cex = 0.8))
p <- p + theme(legend.position="none")
p <- p + theme(text = element_text(size=15))
p
ggsave("sd_sddistance.png", dpi=1000)

